home *** CD-ROM | disk | FTP | other *** search
/ AP Professional Graphics CD-ROM Library / AP Professional Graphics CD-ROM Library.iso / pc / unix / appendix / gemsii / inv_cmap / inv_cmap.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-09-18  |  13.9 KB  |  525 lines
  1. /*
  2.  * This software is copyrighted as noted below.  It may be freely copied,
  3.  * modified, and redistributed, provided that the copyright notice is
  4.  * preserved on all copies.
  5.  *
  6.  * There is no warranty or other guarantee of fitness for this software,
  7.  * it is provided solely "as is".  Bug reports or fixes may be sent
  8.  * to the author, who may or may not act on them as he desires.
  9.  *
  10.  * You may not include this software in a program or other software product
  11.  * without supplying the source, or without informing the end-user that the
  12.  * source is available for no extra charge.
  13.  *
  14.  * If you modify this software, you should include a notice giving the
  15.  * name of the person performing the modification, the date of modification,
  16.  * and the reason for such modification.
  17.  */
  18. /*
  19.  * inv_cmap.c - Compute an inverse colormap.
  20.  *
  21.  * Author:    Spencer W. Thomas
  22.  *         EECS Dept.
  23.  *         University of Michigan
  24.  * Date:    Thu Sep 20 1990
  25.  * Copyright (c) 1990, University of Michigan
  26.  *
  27.  * $Id: inv_cmap.c,v 3.0.1.3 1992/04/30 14:07:28 spencer Exp $
  28.  */
  29.  
  30. #include <math.h>
  31. #include <stdio.h>
  32.  
  33.  
  34. static int bcenter, gcenter, rcenter;
  35. static long gdist, rdist, cdist;
  36. static long cbinc, cginc, crinc;
  37. static unsigned long *gdp, *rdp, *cdp;
  38. static unsigned char *grgbp, *rrgbp, *crgbp;
  39. static gstride, rstride;
  40. static long x, xsqr, colormax;
  41. static int cindex;
  42.  
  43. #ifdef USE_PROTOTYPES
  44. static void maxfill( unsigned long *, long );
  45. static int redloop( void );
  46. static int greenloop( int );
  47. static int blueloop( int );
  48. #else
  49. static void maxfill();
  50. static int redloop();
  51. static int greenloop();
  52. static int blueloop();
  53. #endif
  54.  
  55.  
  56. /*****************************************************************
  57.  * TAG( inv_cmap )
  58.  *
  59.  * Compute an inverse colormap efficiently.
  60.  * Inputs:
  61.  *     colors:        Number of colors in the forward colormap.
  62.  *     colormap:    The forward colormap.
  63.  *     bits:        Number of quantization bits.  The inverse
  64.  *             colormap will have (2^bits)^3 entries.
  65.  *     dist_buf:    An array of (2^bits)^3 long integers to be
  66.  *             used as scratch space.
  67.  * Outputs:
  68.  *     rgbmap:        The output inverse colormap.  The entry
  69.  *             rgbmap[(r<<(2*bits)) + (g<<bits) + b]
  70.  *             is the colormap entry that is closest to the
  71.  *             (quantized) color (r,g,b).
  72.  * Assumptions:
  73.  *     Quantization is performed by right shift (low order bits are
  74.  *     truncated).  Thus, the distance to a quantized color is
  75.  *     actually measured to the color at the center of the cell
  76.  *     (i.e., to r+.5, g+.5, b+.5, if (r,g,b) is a quantized color).
  77.  * Algorithm:
  78.  *     Uses a "distance buffer" algorithm:
  79.  *     The distance from each representative in the forward color map
  80.  *     to each point in the rgb space is computed.  If it is less
  81.  *     than the distance currently stored in dist_buf, then the
  82.  *     corresponding entry in rgbmap is replaced with the current
  83.  *     representative (and the dist_buf entry is replaced with the
  84.  *     new distance).
  85.  *
  86.  *     The distance computation uses an efficient incremental formulation.
  87.  *
  88.  *     Distances are computed "outward" from each color.  If the
  89.  *     colors are evenly distributed in color space, the expected
  90.  *     number of cells visited for color I is N^3/I.
  91.  *     Thus, the complexity of the algorithm is O(log(K) N^3),
  92.  *     where K = colors, and N = 2^bits.
  93.  */
  94.  
  95. /*
  96.  * Here's the idea:  scan from the "center" of each cell "out"
  97.  * until we hit the "edge" of the cell -- that is, the point
  98.  * at which some other color is closer -- and stop.  In 1-D,
  99.  * this is simple:
  100.  *     for i := here to max do
  101.  *         if closer then buffer[i] = this color
  102.  *         else break
  103.  *     repeat above loop with i := here-1 to min by -1
  104.  *
  105.  * In 2-D, it's trickier, because along a "scan-line", the
  106.  * region might start "after" the "center" point.  A picture
  107.  * might clarify:
  108.  *         |    ...
  109.  *               | ...    .
  110.  *              ...        .
  111.  *           ... |      .
  112.  *          .    +         .
  113.  *           .          .
  114.  *            .         .
  115.  *             .........
  116.  *
  117.  * The + marks the "center" of the above region.  On the top 2
  118.  * lines, the region "begins" to the right of the "center".
  119.  *
  120.  * Thus, we need a loop like this:
  121.  *     detect := false
  122.  *     for i := here to max do
  123.  *         if closer then
  124.  *             buffer[..., i] := this color
  125.  *             if !detect then
  126.  *                 here = i
  127.  *                 detect = true
  128.  *         else
  129.  *             if detect then
  130.  *                 break
  131.  *                 
  132.  * Repeat the above loop with i := here-1 to min by -1.  Note that
  133.  * the "detect" value should not be reinitialized.  If it was
  134.  * "true", and center is not inside the cell, then none of the
  135.  * cell lies to the left and this loop should exit
  136.  * immediately.
  137.  *
  138.  * The outer loops are similar, except that the "closer" test
  139.  * is replaced by a call to the "next in" loop; its "detect"
  140.  * value serves as the test.  (No assignment to the buffer is
  141.  * done, either.)
  142.  *
  143.  * Each time an outer loop starts, the "here", "min", and
  144.  * "max" values of the next inner loop should be
  145.  * re-initialized to the center of the cell, 0, and cube size,
  146.  * respectively.  Otherwise, these values will carry over from
  147.  * one "call" to the inner loop to the next.  This tracks the
  148.  * edges of the cell and minimizes the number of
  149.  * "unproductive" comparisons that must be made.
  150.  *
  151.  * Finally, the inner-most loop can have the "if !detect"
  152.  * optimized out of it by splitting it into two loops: one
  153.  * that finds the first color value on the scan line that is
  154.  * in this cell, and a second that fills the cell until
  155.  * another one is closer:
  156.  *      if !detect then        {needed for "down" loop}
  157.  *         for i := here to max do
  158.  *         if closer then
  159.  *             buffer[..., i] := this color
  160.  *             detect := true
  161.  *             break
  162.  *    for i := i+1 to max do
  163.  *        if closer then
  164.  *             buffer[..., i] := this color
  165.  *         else
  166.  *             break
  167.  *
  168.  * In this implementation, each level will require the
  169.  * following variables.  Variables labelled (l) are local to each
  170.  * procedure.  The ? should be replaced with r, g, or b:
  171.  *      cdist:            The distance at the starting point.
  172.  *     ?center:    The value of this component of the color
  173.  *      c?inc:            The initial increment at the ?center position.
  174.  *     ?stride:    The amount to add to the buffer
  175.  *             pointers (dp and rgbp) to get to the
  176.  *             "next row".
  177.  *     min(l):        The "low edge" of the cell, init to 0
  178.  *     max(l):        The "high edge" of the cell, init to
  179.  *             colormax-1
  180.  *     detect(l):        True if this row has changed some
  181.  *                     buffer entries.
  182.  *      i(l):             The index for this row.
  183.  *      ?xx:            The accumulated increment value.
  184.  *      
  185.  *      here(l):        The starting index for this color.  The
  186.  *                      following variables are associated with here,
  187.  *                      in the sense that they must be updated if here
  188.  *                      is changed.
  189.  *      ?dist:            The current distance for this level.  The
  190.  *                      value of dist from the previous level (g or r,
  191.  *                      for level b or g) initializes dist on this
  192.  *                      level.  Thus gdist is associated with here(b)).
  193.  *      ?inc:            The initial increment for the row.
  194.  *      ?dp:            Pointer into the distance buffer.  The value
  195.  *                      from the previous level initializes this level.
  196.  *      ?rgbp:            Pointer into the rgb buffer.  The value
  197.  *                      from the previous level initializes this level.
  198.  * 
  199.  * The blue and green levels modify 'here-associated' variables (dp,
  200.  * rgbp, dist) on the green and red levels, respectively, when here is
  201.  * changed.
  202.  */
  203.  
  204. void
  205. inv_cmap( colors, colormap, bits, dist_buf, rgbmap )
  206. int colors, bits;
  207. unsigned char *colormap[3], *rgbmap;
  208. unsigned long *dist_buf;
  209. {
  210.     int nbits = 8 - bits;
  211.  
  212.     colormax = 1 << bits;
  213.     x = 1 << nbits;
  214.     xsqr = 1 << (2 * nbits);
  215.  
  216.     /* Compute "strides" for accessing the arrays. */
  217.     gstride = colormax;
  218.     rstride = colormax * colormax;
  219.  
  220.     maxfill( dist_buf, colormax );
  221.  
  222.     for ( cindex = 0; cindex < colors; cindex++ )
  223.     {
  224.     /*
  225.      * Distance formula is
  226.      * (red - map[0])^2 + (green - map[1])^2 + (blue - map[2])^2
  227.      *
  228.      * Because of quantization, we will measure from the center of
  229.      * each quantized "cube", so blue distance is
  230.      *     (blue + x/2 - map[2])^2,
  231.      * where x = 2^(8 - bits).
  232.      * The step size is x, so the blue increment is
  233.      *     2*x*blue - 2*x*map[2] + 2*x^2
  234.      *
  235.      * Now, b in the code below is actually blue/x, so our
  236.      * increment will be 2*(b*x^2 + x^2 - x*map[2]).  For
  237.      * efficiency, we will maintain this quantity in a separate variable
  238.      * that will be updated incrementally by adding 2*x^2 each time.
  239.      */
  240.     /* The initial position is the cell containing the colormap
  241.      * entry.  We get this by quantizing the colormap values.
  242.      */
  243.     rcenter = colormap[0][cindex] >> nbits;
  244.     gcenter = colormap[1][cindex] >> nbits;
  245.     bcenter = colormap[2][cindex] >> nbits;
  246.  
  247.     rdist = colormap[0][cindex] - (rcenter * x + x/2);
  248.     gdist = colormap[1][cindex] - (gcenter * x + x/2);
  249.     cdist = colormap[2][cindex] - (bcenter * x + x/2);
  250.     cdist = rdist*rdist + gdist*gdist + cdist*cdist;
  251.  
  252.     crinc = 2 * ((rcenter + 1) * xsqr - (colormap[0][cindex] * x));
  253.     cginc = 2 * ((gcenter + 1) * xsqr - (colormap[1][cindex] * x));
  254.     cbinc = 2 * ((bcenter + 1) * xsqr - (colormap[2][cindex] * x));
  255.  
  256.     /* Array starting points. */
  257.     cdp = dist_buf + rcenter * rstride + gcenter * gstride + bcenter;
  258.     crgbp = rgbmap + rcenter * rstride + gcenter * gstride + bcenter;
  259.  
  260.     (void)redloop();
  261.     }
  262. }
  263.  
  264. /* redloop -- loop up and down from red center. */
  265. static int
  266. redloop()
  267. {
  268.     int detect;
  269.     int r;
  270.     int first;
  271.     long txsqr = xsqr + xsqr;
  272.     static long rxx;
  273.  
  274.     detect = 0;
  275.  
  276.     /* Basic loop up. */
  277.     for ( r = rcenter, rdist = cdist, rxx = crinc,
  278.       rdp = cdp, rrgbp = crgbp, first = 1;
  279.       r < colormax;
  280.       r++, rdp += rstride, rrgbp += rstride,
  281.       rdist += rxx, rxx += txsqr, first = 0 )
  282.     {
  283.     if ( greenloop( first ) )
  284.         detect = 1;
  285.     else if ( detect )
  286.         break;
  287.     }
  288.     
  289.     /* Basic loop down. */
  290.     for ( r = rcenter - 1, rxx = crinc - txsqr, rdist = cdist - rxx,
  291.       rdp = cdp - rstride, rrgbp = crgbp - rstride, first = 1;
  292.       r >= 0;
  293.       r--, rdp -= rstride, rrgbp -= rstride,
  294.       rxx -= txsqr, rdist -= rxx, first = 0 )
  295.     {
  296.     if ( greenloop( first ) )
  297.         detect = 1;
  298.     else if ( detect )
  299.         break;
  300.     }
  301.     
  302.     return detect;
  303. }
  304.  
  305. /* greenloop -- loop up and down from green center. */
  306. static int
  307. greenloop( restart )
  308. int restart;
  309. {
  310.     int detect;
  311.     int g;
  312.     int first;
  313.     long txsqr = xsqr + xsqr;
  314.     static int here, min, max;
  315.     static long ginc, gxx, gcdist;    /* "gc" variables maintain correct */
  316.     static unsigned long *gcdp;        /*  values for bcenter position, */
  317.     static unsigned char *gcrgbp;    /*  despite modifications by blueloop */
  318.                     /*  to gdist, gdp, grgbp. */
  319.  
  320.     if ( restart )
  321.     {
  322.     here = gcenter;
  323.     min = 0;
  324.     max = colormax - 1;
  325.     ginc = cginc;
  326.     }
  327.  
  328.     detect = 0;
  329.  
  330.     /* Basic loop up. */
  331.     for ( g = here, gcdist = gdist = rdist, gxx = ginc,
  332.       gcdp = gdp = rdp, gcrgbp = grgbp = rrgbp, first = 1;
  333.       g <= max;
  334.       g++, gdp += gstride, gcdp += gstride, grgbp += gstride, gcrgbp += gstride,
  335.       gdist += gxx, gcdist += gxx, gxx += txsqr, first = 0 )
  336.     {
  337.     if ( blueloop( first ) )
  338.     {
  339.         if ( !detect )
  340.         {
  341.         /* Remember here and associated data! */
  342.         if ( g > here )
  343.         {
  344.             here = g;
  345.             rdp = gcdp;
  346.             rrgbp = gcrgbp;
  347.             rdist = gcdist;
  348.             ginc = gxx;
  349.         }
  350.         detect = 1;
  351.         }
  352.     }
  353.     else if ( detect )
  354.     {
  355.         break;
  356.     }
  357.     }
  358.     
  359.     /* Basic loop down. */
  360.     for ( g = here - 1, gxx = ginc - txsqr, gcdist = gdist = rdist - gxx,
  361.       gcdp = gdp = rdp - gstride, gcrgbp = grgbp = rrgbp - gstride,
  362.       first = 1;
  363.       g >= min;
  364.       g--, gdp -= gstride, gcdp -= gstride, grgbp -= gstride, gcrgbp -= gstride,
  365.       gxx -= txsqr, gdist -= gxx, gcdist -= gxx, first = 0 )
  366.     {
  367.     if ( blueloop( first ) )
  368.     {
  369.         if ( !detect )
  370.         {
  371.         /* Remember here! */
  372.         here = g;
  373.         rdp = gcdp;
  374.         rrgbp = gcrgbp;
  375.         rdist = gcdist;
  376.         ginc = gxx;
  377.         detect = 1;
  378.         }
  379.     }
  380.     else if ( detect )
  381.     {
  382.         break;
  383.     }
  384.     }
  385.     
  386.  
  387.     return detect;
  388. }
  389.  
  390. /* blueloop -- loop up and down from blue center. */
  391. static int
  392. blueloop( restart )
  393. int restart;
  394. {
  395.     int detect;
  396.     register unsigned long *dp;
  397.     register unsigned char *rgbp;
  398.     register long bdist, bxx;
  399.     register int b, i = cindex;
  400.     register long txsqr = xsqr + xsqr;
  401.     register int lim;
  402.     static int here, min, max;
  403.     static long binc;
  404.  
  405.     if ( restart )
  406.     {
  407.     here = bcenter;
  408.     min = 0;
  409.     max = colormax - 1;
  410.     binc = cbinc;
  411.     }
  412.  
  413.     detect = 0;
  414.  
  415.     /* Basic loop up. */
  416.     /* First loop just finds first applicable cell. */
  417.     for ( b = here, bdist = gdist, bxx = binc, dp = gdp, rgbp = grgbp, lim = max;
  418.       b <= lim;
  419.       b++, dp++, rgbp++,
  420.       bdist += bxx, bxx += txsqr )
  421.     {
  422.     if ( *dp > bdist )
  423.     {
  424.         /* Remember new 'here' and associated data! */
  425.         if ( b > here )
  426.         {
  427.         here = b;
  428.         gdp = dp;
  429.         grgbp = rgbp;
  430.         gdist = bdist;
  431.         binc = bxx;
  432.         }
  433.         detect = 1;
  434.         break;
  435.     }
  436.     }
  437.     /* Second loop fills in a run of closer cells. */
  438.     for ( ;
  439.       b <= lim;
  440.       b++, dp++, rgbp++,
  441.       bdist += bxx, bxx += txsqr )
  442.     {
  443.     if ( *dp > bdist )
  444.     {
  445.         *dp = bdist;
  446.         *rgbp = i;
  447.     }
  448.     else
  449.     {
  450.         break;
  451.     }
  452.     }
  453.     
  454.     /* Basic loop down. */
  455.     /* Do initializations here, since the 'find' loop might not get
  456.      * executed. 
  457.      */
  458.     lim = min;
  459.     b = here - 1;
  460.     bxx = binc - txsqr;
  461.     bdist = gdist - bxx;
  462.     dp = gdp - 1;
  463.     rgbp = grgbp - 1;
  464.     /* The 'find' loop is executed only if we didn't already find
  465.      * something.
  466.      */
  467.     if ( !detect )
  468.     for ( ;
  469.           b >= lim;
  470.           b--, dp--, rgbp--,
  471.           bxx -= txsqr, bdist -= bxx )
  472.     {
  473.         if ( *dp > bdist )
  474.         {
  475.         /* Remember here! */
  476.         /* No test for b against here necessary because b <
  477.          * here by definition.
  478.          */
  479.         here = b;
  480.         gdp = dp;
  481.         grgbp = rgbp;
  482.         gdist = bdist;
  483.         binc = bxx;
  484.         detect = 1;
  485.         break;
  486.         }
  487.     }
  488.     /* The 'update' loop. */
  489.     for ( ;
  490.       b >= lim;
  491.       b--, dp--, rgbp--,
  492.       bxx -= txsqr, bdist -= bxx )
  493.     {
  494.     if ( *dp > bdist )
  495.     {
  496.         *dp = bdist;
  497.         *rgbp = i;
  498.     }
  499.     else
  500.     {
  501.         break;
  502.     }
  503.     }
  504.  
  505.  
  506.     /* If we saw something, update the edge trackers. */
  507.  
  508.     return detect;
  509. }
  510.  
  511. static void
  512. maxfill( buffer, side )
  513. unsigned long *buffer;
  514. long side;
  515. {
  516.     register unsigned long maxv = ~0L;
  517.     register long i;
  518.     register unsigned long *bp;
  519.  
  520.     for ( i = side * side * side, bp = buffer;
  521.       i > 0;
  522.       i--, bp++ )
  523.     *bp = maxv;
  524. }
  525.